/*
 * nr_rk4_phase_lattice.cpp
 *
 * Copyright 2009-2012 Karsten Ahnert
 * Copyright 2009-2012 Mario Mulansky
 *
 * Distributed under the Boost Software License, Version 1.0.
 * (See accompanying file LICENSE_1_0.txt or
 * copy at http://www.boost.org/LICENSE_1_0.txt)
 */

#include <boost/array.hpp>

#include "rk_performance_test_case.hpp"

#include "phase_lattice.hpp"

const size_t dim = 1024;

typedef boost::array<double, dim> state_type;

template <class System, typename T, size_t dim>
void rk4_step(const System sys, boost::array<T, dim> &x, const double t,
              const double dt) {  // fast rk4 implementation adapted from the
                                  // book 'Numerical Recipes'
  size_t i;
  const double hh = dt * 0.5;
  const double h6 = dt / 6.0;
  const double th = t + hh;
  boost::array<T, dim> dydx, dym, dyt, yt;

  sys(x, dydx, t);

  for (i = 0; i < dim; i++) yt[i] = x[i] + hh * dydx[i];

  sys(yt, dyt, th);
  for (i = 0; i < dim; i++) yt[i] = x[i] + hh * dyt[i];

  sys(yt, dym, th);
  for (i = 0; i < dim; i++) {
    yt[i] = x[i] + dt * dym[i];
    dym[i] += dyt[i];
  }
  sys(yt, dyt, t + dt);
  for (i = 0; i < dim; i++) x[i] += h6 * (dydx[i] + dyt[i] + 2.0 * dym[i]);
}

class nr_wrapper {
 public:
  void reset_init_cond() {
    for (size_t i = 0; i < dim; ++i)
      m_x[i] = 2.0 * 3.1415927 * rand() / RAND_MAX;
    m_t = 0.0;
  }

  inline void do_step(const double dt) {
    rk4_step(phase_lattice<dim>(), m_x, m_t, dt);
  }

  double state(const size_t i) const { return m_x[i]; }

 private:
  state_type m_x;
  double m_t;
};

int main() {
  srand(12312354);

  nr_wrapper stepper;

  run(stepper, 10000, 1E-6);
}
